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We use molecular dynamics simulations to study a model of the gelation transition with a dynamic 
bond forming procedure. After establishing evidence for 3D percolation as the static universality 
class, we turn our attention to the dynamics of clusters at the gelation point, focusing in particular 
on the behavior of the diffusion constant D(s) as a function of cluster size s. We find a very clear 
■ power law behavior D(s) ~ s~ k with a non-trivial exponent k « 0.69. 
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PACS numbers: 82.70.Gg, 61.43.Hv, 36.40.Sx 



^ . I. INTRODUCTION 

O 



Gelation is a phenomenon that continues to receive wide interest in the physics community (l], |^, || . It is a problem 
with relevance both from a fundamental physics point of view, but also from a more application oriented view due 
to the many uses of gels in industry. Despite a long-standing effort the underlying physics has still not received a 
complete explanation, in particular as far as the dynamical quantities are concerned. 

In a solution of polyfunctional monomers irreversible chemical bonds may be formed via some external influence, 
e.g. irradiation with 7-rays |5j. As the density p of bonds increases, the resulting macromolecules grow larger 
and larger, until at a critical density p c (the gelation point) there is a cluster spanning the entire system, and a gel 
has formed. There are a number of static/thermodynamic quantities characterizing this transition, e.g. cluster size 
[ distribution, largest cluster, correlation length and radius of gyration, that in many cases are believed to be well 
described by the 3 dimensional (3-D) percolation universality class ||. Besides the appearance of a wide distribution 
of cluster sizes, the presence of the gel fundamentally changes many of the dynamic properties as well. For example 
CD . it has been found that the viscosity rj diverges as (p c — p)~ s when p — > p c , but there is still no general agreement on 
the value of the exponent s. However there are many other quantities that probe the unique dynamics at the gelation 
-H , point, and here we consider, among other things, the diffusion of clusters of different sizes. Diffusion in general is well 
known to be sensitive to geometry and topology as well as to the interactions with the environment, and in the case 
of gelation we expect the motion of the individual macromolecules to be strongly affected by their interaction with 
other macromolecules with a wide distribution of sizes. 

Here we consider a model for polycondensation: using molecular dynamics we simulate a solution of hexafunctional 
monomers interacting via the Lennard- Jones potential and allow the particles to form permanent chemical bonds if 
{SJ ' they get close enough together. After defining the model in Sec. II, we determine the static universality class by 
considering the behavior of several characteristic structural quantities aided by finite size scaling in Sec. [ffl]. Here we 
' also measure the radius of gyration of the clusters as a function of cluster size at the gelation point allowing us to 
determine the fractal dimension. In Sec. |fv| we focus on the behavior of the self-diffusion constant of clusters at the 



gelation point, a quantity which has hitherto received little attention (see however 0). Finally in Sec. [v]we give our 
conclusions. 

II. MODEL 

Our system is composed of N = L 3 (L — 10, 15 and 20) particles interacting pairwise through the shifted Lennard- 
Jones potential 

u{r)= iu LJ (r)-U LJ (2.5a), r < 2.5a 
1 otherwise, 

where ULj(r) = 4e((cr/r) 12 — (a/r) 6 ). All of our simulations are 3D constant energy simulations corresponding to 
an average temperature of ksT/e « 1 and density <& = 0.8<r~ 3 . These choices ensure that the system is in the 
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FIG. 1: Fraction of systems W(L,p) percolating in the £- direct ion as a function of p and for 3 different system sizes as 
indicated on the plot. The lines are guides for the eye. We estimate p c = 0.2565. 

liquid-phase region of the phase-diagram |H. We use periodic boundary conditions and a time step of magnitude 
dt = 0. 004 -y/TOer 2 Je. From a typical equilibrium state of this liquid we let the particles form permanent chemical 
bonds if they come closer than r c — 2 1 / 6 cr w 1.12 (coinciding with the minimum of U{r)), and the corresponding 
bond interaction is represented by a harmonic oscillator potential Uh aim (r) = l/2kr 2 : in our simulations we take 
ka 2 /e = 120. Note that this way of adding bonds breaks energy conservation; indeed we actually pump energy into 
the system when adding bonds. With this bonding procedure cross linking is very fast — the average distance between 
the particles is comparable to r c , so a large number of particles will be available for bonding at a given instant. Each 
particle can bond to a total of / = 6 other particles, and the cross link density p is then given in terms of the number 
of bonds n as p = 2n/ fN. Any number of particles, if fulfilling the conditions above, can be cross linked per time 
step, but we halt the bond formation when p reaches a predetermined value. 

III. STATIC UNIVERSALITY CLASS 

We begin by examining the geometric structure resulting from the cross linking procedure described above. First 
we need to know the location of the geometric percolation point p c , and to this end we measure the fraction W(L,p) 
of percolating (in one direction) systems of size N — L 3 at a given cross link density p. In this context a system 
is considered percolating when a cluster connects a particle in the central computational box with its image in an 
adjacent box, and the corresponding cluster is called a spanning cluster. In the limit N — > oo, W(L,p) becomes 
the Heaviside step function 6{p — p c ). However at finite N for any p > below p Cl some realizations will contain a 
spanning cluster, and thus W(L,p) > 0. Similar reasoning allows us to conclude that W(L,p) < 1 for p > p c . If 
we make the plausible assumption that W(L,p) converges monotonically to 9(p — p c ) for fixed p, it follows that for 
L\ > L 2 , W{L\,p) < W{L 2l p) if p < Pc, and likewise W{L\,p) > W(L 2 ,p) for any p > p c . Therefore at and only at 
p c the curves of W(L,p) coincide, i.e. W(L,p c ) = const. In Fig. |l| we plot W{L,p) obtained from 5000 realizations as 
a function of p for L = 10, 15 and 20, and we see that all the curves cross approximately at the same point: the value 
of p c determined this way is p c = 0.2565. 

Finite size scaling theory predicts that W(L,p) does not depend on L and p separately but only on the combination 
£/£ (and the sign of p — p c ) where £ = \p — p c \~ u is the correlation length and v the correlation length exponent ||. 
Thus we may write 

W(L lP ) = f(L 1 ^(p-p c )), (2) 
where f{x) is a scaling function with the following limiting values: 



lim f(x) = 1 and lim f(x) = 0. 



(3) 
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FIG. 3: Here we plot the cluster size distribution n(s) as a function of s for the same 3 system sizes at p = p c and based on 
5000 samples. The straight line is a fit to a power law to the data from the L — 20 system. 



To test this hypothesis we replot the data for W(L,p) from Fig. [j] in Fig. |] as a function of L}/ V {p — p c ) with 
p c = 0.2565 as determined above and v = 0.9. The curves collapse quite convincingly onto the same master curve 
with the limiting values given in Eq. (|3|), and the results are thus consistent with finite scaling theory. Moreover, the 
value of v is in very good agreement with the correlation length exponent 0.88 of 3D percolation theory [|| g]. Thus 
we believe that in spite of using an off-lattice dynamic bonding method from an equilibrium liquid state, this systems 
belongs to the universality class of 3D percolation. To further substantiate that claim, we have also considered the 
cluster size distribution at the critical point p = p c in Fig. ^. According to standard percolation theory the number 
n(s) of clusters of size s is a power law n(s) ~ s _T with an exponent r « 2.18 jlfj[| . Fitting the data from the largest 
system size to a power law we find the exponent r w 2.19, which is within 0.5% of the 3D percolation value quoted 
above. 

For the universality class of 3D percolation there are in principle only two independent exponents, but to establish 
even more confidence in our conclusions and because we need part of the data later, we determine two more exponents. 
The weight average molecular mass M w for p < p c defines the exponent 7 by the relation M w ~ (p c —p) ' . Again, 
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FIG. 4: Scaling plot of the weight average molecular weight M w . The quality of the data collapse confirms 7 = 1.74 in 
accordance with the 3D percolation value. 



for finite systems, finite size scaling theory predicts the scaling form 

M w =L^g(L 1 ^(p c -p)), (4) 
where g(x) is a scaling function with the following asymptotic behavior: 

g(x) ~ X °° (5) 

I const, x — > 0. 

Therefore we compute M ffi as a function of p for different system sizes, and in Fig. [| we plot the results in the form 
M w jL 1 l v versus L x l v (j) c — p) with 7 = 1.74 being the expected 3D percolation value and v and p c as determined 
previously. Again there is a very nice data collapse reinforcing our statement about the validity of 3D percolation 
exponents. 

The last static exponent we wish to consider is the fractal dimension D of the clusters at the gelation point. During 
the simulations we maintain a list of all particles belonging to a given cluster, and this enables us to calculate the 
radius of gyration of a given cluster. At a given instant R g (s) (s > 2) is found from |ll| 

V ' i,J = l 

where the sum is over all particles in the cluster and rj is the position of particle i at that instant. To calculate its 
thermal average we average over time. However we have also averaged over 100 samples to get better statistics for 
the individual cluster sizes. We find that the results are independent of the system size (up to a cut-off) and do not 
significantly depend on the value of p close to p c , and therefore we plot here only the results from the largest system, 
i.e. N = 20 3 . We focus on the properties of the system at the percolation point p = p c and the corresponding data 
are plotted in Fig. g. 

There is a clear power law behavior Rg(s) ~ s 2 l D , and from a least squares fit we deduce the value D = 2.38. Here 
D is the fractal dimension of the clusters (since s ~ R®) at the gelation point, and for this exponent 3D percolation 
theory predicts the value D « 2.5. Our value is about 5% lower, but this small discrepancy may be due to lack of 
data for sufficiently large s. However it may also be due to the interaction of the polymers with one another and with 
the left over monomers. 

Therefore regarding the static exponents we presented evidence that our model of the gelation transition, which 
uses a Lennard- Jones force-shifted intermolecular potential and a liquid state bond forming procedure, is consistent 
with the universality class of 3D percolation. 
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FIG. 5: Radius of gyrations R g (squared) as a function of cluster size s for L — 20 and p = p c . 




Timcstcps Timcstcps 

(a) (b) 

FIG. 6: Mean square displacements divided by time for the individual clusters of the N = 20 3 system at p = p c . (a) Clusters 
of size 1 ... 29 from top to bottom. The asymptotic convergence to horizontal lines indicates normal diffusion, c.f. Eq. (Q). (b) 
Close up of (a) for the 4 largest clusters. They have all reached an almost constant value, except perhaps the largest cluster, 
which therefore represents the "worst case" . 



IV. CLUSTER DIFFUSION 



We now turn to the measurement of a dynamic cluster specific quantity, namely the diffusion coefficient D(s) as a 
function of cluster size s at the percolation threshold p — p c . We measure the mean square fluctuations of the center 
of mass position R^ M (t) for each cluster of size s, and for large times we find a linear behavior 

((R™(t)-R™(0)f)~6D(s)t, (7) 

allowing us to extract the diffusion constant D(s) as a function of cluster mass s. Again the average (■ • • } is performed 
over both time (thermal average) and over realizations, and here we have found it sufficient to use 80 samples. 

To demonstrate the asymptotic normal diffusion we plot for the N = 20 3 system in Fig. ||.a {(R™ (t) - R^ M (0)) 2 ) /t 
as a function of time steps for several cluster sizes s — 1 ... 29, and part (b) shows a magnification for the largest 
clusters s = 26 . . . 29. As can be seen all the curves approach constant values for large times, and it is from these 
values that we calculate the diffusion constant using Eq. (Q). We have checked that the results below do not change 
if we extend the length of the sampling time to much longer times (5000 timesteps). 

For monodisperse polymers in a dilute solution, Rouse dynamics predicts the power law dependence D(s) ~ s . 
Taking hydrodynamic interactions into account, Zimm dynamics predicts D(s) ~ R~ x ~ s~ x / D = s~ 0A p"2] ]. For the 
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FIG. 7: Diffusion constant D(s) as a function of the mass s (or size) of the cluster for the N = 20 3 system at p = p c . The 
solid line is a power law s -0 ' 69 obtained from a least squares fit to the data in the region [3: 30]. 
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FIG. 



Diffusion constant as a function of radius of gyration for N = 20 3 and p = p c 



present model we observe a very clear power law D(s) ~ s~ k as can be seen in Fig. fj], where we have plotted D(s) as 
a function of s. We consider again only the largest system size N = 20 3 at p = p c , and D(s) is found from Eq. (]?]) 
based on the data from Fig. |6| The exponent of the power law is approximately k « 0.69, and is thus described 
neither by pure Rouse nor Zimm dynamics. Instead the exponent is intermediate between the values for Rouse and 
Zimm dynamics. 

Having measured the radius of gyration in Fig. |5|, we can convert the result above into a dependence of D on 
R g : D(Rg) ~ s~ - 69 ^ R~ o mD = R g M , and the corresponding data are plotted in Fig. [|. Also shown is a power 
law fit to the data yielding an exponent value of 1.63. Again we emphasize the clear deviation from a Stokes' type 
law, D{Rg) r~j R^ 1 , as obtained from Zimm dynamics for polymers. Moreover our result is in sharp contrast to a 
corresponding result in a recent paper [jl3| for a lattice model of gelation in which it was found that D{R g ) ~ R~ 2A 
using bond fluctuation dynamics. 
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V. CONCLUSIONS 



Using extensive molecular dynamics simulations of a system of Lennard- Jones particles we studied a specific model 
for the gelation transition. An important feature of our model is that particles are allowed to form permanent chemical 
bonds from the Lennard- Jones liquid state when they get close together. For this model we determined a number of 
exponents characterizing the geometry and polydispersity of the system, and this information allowed us to conclude 
that 3D percolation theory is the appropriate static universality class. We went on to discuss the radius of gyration 
of the clusters at the gel point resulting in the determination of their fractal dimension. Afterwards we focused on 
the self-diffusion constant of individual clusters as a function of their size. We found strong indications for a power 
law behavior with an exponent described neither by Rouse nor Zimm-dynamics. In future work we intend to study 
e.g. the shear viscosity of this model in order to measure the corresponding dynamic exponents. As mentioned in the 
introduction there is considerable disagreement as to the values of these exponents, and numerical results on models 
such as the present one, will allow us to discriminate between some of the proposed dynamic universality classes for 
gelation. 
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